Result_OriginalData

Library

library(tidyverse)
library(here)      
library(readxl)    
library(janitor)   
library(stringr)   
library(tidyr)     
library(dplyr)
library(knitr)
library(ggplot2)
library(factoextra)
library(gt)
library(naniar)
library(ggcorrplot)
library(corrplot)
library(GGally)
library(randomForest)
library(e1071)  
library(caret)
library(plotly)
library(stargazer)
library(broom)
library(kableExtra)
library(gridExtra)
library(patchwork)
library(minerva)
library(mgcv)
library(MASS)
library(glmnet)
library(purrr)
library(htmltools)
library(insight)

Data Import

raw_bp <- read.csv(here("ZIPallFiles","AHSnpa11bp.csv"), header=TRUE)
raw_bb <- read.csv(here("ZIPallFiles","AHSnpa11bb.csv"), header=TRUE)
raw_bsp <- read.csv(here("ZIPallFiles","inp13bsp.csv"), header=TRUE)
raw_bcn <- read.csv(here("ZIPallFiles","inp13bcn.csv"), header=TRUE)
raw_bhh <- read.csv(here("ZIPallFiles","inp13bhh.csv"), header=TRUE)

Rename and Created

bp_bb <- inner_join(raw_bp, raw_bb, by = "ABSPID") %>%
  rename(HypertensionStatus = HYPBC,
         BodyMass = SABDYMS,
         SocialStatus = SF2SA1QN,
         PersonID = ABSPID,
         Age = AGEC,
         Gender = SEX,
         BMI = BMISC,
         Systolic = SYSTOL,
         Diastolic = DIASTOL,
         SmokeStatus = SMKSTAT,
         AlcoholPercentage_Day1 = ALCPER1,
         AlcoholPercentage_Day2 = ALCPER2,
         Potassium_Day1 = POTAST1,
         Potassium_Day2 = POTAST2,
         Sodium_Day1 = SODIUMT1,
         Sodium_Day2 = SODIUMT2,
         FibreEnergy_Day1 = FIBRPER1,
         FibreEnergy_Day2 = FIBRPER2,
         CarbohydrateEnergy_Day1 = CHOPER1,
         CarbohydrateEnergy_Day2 = CHOPER2,
         EnergyBMR_Day1 = EIBMR1,
         EnergyBMR_Day2 = EIBMR2) %>%
   mutate(Identity = factor(0))

bsp_bhh<- left_join(raw_bsp, raw_bhh,by = c("ABSHID")) %>%
  rename(BodyMass = SABDYMS,
         SocialStatus = SF2SA1DB,
         HouseID = ABSHID,
         Age = AGEEC,
         Gender = SEX,
         BMI = BMISC,
         Systolic = SYSTOL,
         Diastolic = DIASTOL,
         SmokeStatus = SMKSTAT,
         AlcoholPercentage_Day1 = ALCPER1,
         AlcoholPercentage_Day2 = ALCPER2,
         Potassium_Day1 = POTAST1,
         Potassium_Day2 = POTAST2,
         Sodium_Day1 = SODIUMT1,
         Sodium_Day2 = SODIUMT2,
         FibreEnergy_Day1 = FIBRPER1,
         FibreEnergy_Day2 = FIBRPER2,
         CarbohydrateEnergy_Day1 = CHOPER1,
         CarbohydrateEnergy_Day2 = CHOPER2,
         EnergyBMR_Day1 = EIBMR1,
         EnergyBMR_Day2 = EIBMR2)%>%
   mutate(Identity = factor(1))

bcn <- raw_bcn %>%
  rename(HouseID = ABSHID,
         MedicalCondition = ICD10ME)

First Filter

bp_bb <- bp_bb %>%
  mutate(EnergyBMR_Day1 = na_if(EnergyBMR_Day1, 997) %>% na_if(998),
         EnergyBMR_Day2 = na_if(EnergyBMR_Day2, 997) %>% na_if(998),
         EnergyBMR = (EnergyBMR_Day1 + EnergyBMR_Day2)/2) %>%
  filter(
    Age >= 18, 
    HypertensionStatus == 5,
    BodyMass != 4,
    (is.na(EnergyBMR) | EnergyBMR >= 0.9))

bsp_bhh <- bsp_bhh %>% 
  mutate(EnergyBMR_Day1 = na_if(EnergyBMR_Day1, 997) %>% na_if(998),
         EnergyBMR_Day2 = na_if(EnergyBMR_Day2, 997) %>% na_if(998),
         EnergyBMR = (EnergyBMR_Day1 + EnergyBMR_Day2)/2) %>%
  filter(
    Age >= 18, 
    BodyMass != 4,
    !HouseID %in% bcn$HouseID[bcn$MedicalCondition %in% c(7, 20)],
    (is.na(EnergyBMR) | EnergyBMR >= 0.9))

Merge

raw_data <- bind_rows(bp_bb,bsp_bhh)

Variable Selection

selected_data <- raw_data %>%
  dplyr::select(PersonID, HouseID, SocialStatus, Age, Gender, BMI, Systolic, Diastolic, SmokeStatus, Identity,
         AlcoholPercentage_Day1, AlcoholPercentage_Day2,
         Potassium_Day1, Potassium_Day2, Sodium_Day1, Sodium_Day2,
         FibreEnergy_Day1, FibreEnergy_Day2, 
         CarbohydrateEnergy_Day1, CarbohydrateEnergy_Day2)

Type Check 1

data.frame(Variable = names(selected_data), Type = sapply(selected_data, class)) %>%
  gt() %>%
  tab_header(
    title = "Variable Names and Their Types"
  ) %>%
 tab_caption(caption = md("Table 1: Variable types table"))

Table 1: Variable types table

Variable Names and Their Types
Variable	Type
PersonID	character
HouseID	character
SocialStatus	integer
Age	integer
Gender	integer
BMI	numeric
Systolic	integer
Diastolic	integer
SmokeStatus	integer
Identity	factor
AlcoholPercentage_Day1	numeric
AlcoholPercentage_Day2	numeric
Potassium_Day1	numeric
Potassium_Day2	numeric
Sodium_Day1	numeric
Sodium_Day2	numeric
FibreEnergy_Day1	numeric
FibreEnergy_Day2	numeric
CarbohydrateEnergy_Day1	numeric
CarbohydrateEnergy_Day2	numeric

Na values + New Variables + Type Conversion

selected_data <- selected_data %>%
  mutate(
    across(c(Gender, SocialStatus, SmokeStatus), as.factor),
    BMI = na_if(BMI, 98) %>% na_if(99),
    Systolic = na_if(Systolic, 0) %>% na_if(998) %>% na_if(999),
    Diastolic = na_if(Diastolic, 0) %>% na_if(998) %>% na_if(999),
    Potassium = (Potassium_Day1 + Potassium_Day2) / 2, 
    Sodium = (Sodium_Day1 + Sodium_Day2) / 2,
    FibreEnergy = (FibreEnergy_Day1 + FibreEnergy_Day2)/2,
    CarbohydrateEnergy = (CarbohydrateEnergy_Day1 + CarbohydrateEnergy_Day2) / 2,
    AlcoholPercentage = (AlcoholPercentage_Day1 + AlcoholPercentage_Day2) /2,
    SodiumPotassiumRatio = Sodium / Potassium
    ) %>% 
  dplyr::select(PersonID,HouseID, SocialStatus, Age, Gender, BMI, Systolic, Diastolic, SmokeStatus, AlcoholPercentage, SodiumPotassiumRatio, FibreEnergy, CarbohydrateEnergy, Identity)

Type Check 2

data.frame(Variable = names(selected_data), Type = sapply(selected_data, class)) %>%
  gt() %>%
  tab_header(
    title = "Variable Names and Their Types"
  ) %>%
 tab_caption(caption = md("Table 2: Variable types table after correction"))

Table 2: Variable types table after correction

Variable Names and Their Types
Variable	Type
PersonID	character
HouseID	character
SocialStatus	factor
Age	integer
Gender	factor
BMI	numeric
Systolic	integer
Diastolic	integer
SmokeStatus	factor
AlcoholPercentage	numeric
SodiumPotassiumRatio	numeric
FibreEnergy	numeric
CarbohydrateEnergy	numeric
Identity	factor

Duplicate Values

selected_data %>% 
  dplyr::select(-PersonID) %>% 
  distinct() %>% 
  dim() %>% 
  tibble::tibble(Dimension = c("Rows", "Columns"), Count = .) %>% 
  gt() %>%
  tab_header(title = "Dimensions of Unique Tibble") %>%
  tab_caption(caption = md("Table 3: Dimensions of Unique Data Table"))

Table 3: Dimensions of Unique Data Table

Dimensions of Unique Tibble
Dimension	Count
Rows	8449
Columns	13

Categrocial Distribution

plot_categorical_distribution <- function(data) {
  categorical_vars <- data %>%
    dplyr::select(where(is.factor) | where(is.character)) %>%
    dplyr::select(-PersonID, -HouseID)
  figure_counter <- 1
  for (var in names(categorical_vars)) {
    non_na_data <- data %>%
      filter(!is.na(.data[[var]]))
    if (nrow(non_na_data) == 0) next
    
    figure_title <- paste("Figure 3.", figure_counter, ": Proportion of Categories in", var, sep = "")
    
    p <- ggplot(non_na_data, aes(x = .data[[var]])) +       
      geom_bar(aes(y = (..count..) / sum(..count..), fill = .data[[var]]), color = "black") +       
      scale_fill_brewer(palette = "Set3") +        
      labs(title = figure_title, x = "Category", y = "Proportion") +       
      theme_minimal() +       
      theme(axis.text.x = element_text(angle = 45, hjust = 1))
    print(p)
    figure_counter <- figure_counter + 1
    }
  }

plot_categorical_distribution(selected_data) 

Merge Variable

selected_data <- selected_data %>% 
  mutate(SmokeStatus = recode_factor(SmokeStatus, "2" = "1", "3" = "1","4" = "2","5" = "3"))

Outliers and New Variables

normal_ranges <- list(Systolic = c(90, 180), Diastolic = c(60, 120), BMI = c(16, 47.5), 
                      SodiumPotassiumRatio = c(0, 4), AlcoholPercentage = c(0, 85), 
                      CarbohydrateEnergy = c(0, 100), FibreEnergy = c(0, 100))

selected_data <- selected_data %>%
  filter(
    (is.na(Systolic) | (Systolic >= normal_ranges$Systolic[1] & Systolic <= normal_ranges$Systolic[2])),
    (is.na(Diastolic) | (Diastolic >= normal_ranges$Diastolic[1] & Diastolic <= normal_ranges$Diastolic[2])),
    (is.na(BMI) | (BMI >= normal_ranges$BMI[1] & BMI <= normal_ranges$BMI[2])),
    (is.na(AlcoholPercentage) | (AlcoholPercentage >= normal_ranges$AlcoholPercentage[1] & AlcoholPercentage <= normal_ranges$AlcoholPercentage[2])),
    (is.na(CarbohydrateEnergy) | (CarbohydrateEnergy >= normal_ranges$CarbohydrateEnergy[1] & CarbohydrateEnergy <= normal_ranges$CarbohydrateEnergy[2])),
    (is.na(FibreEnergy) | (FibreEnergy >= normal_ranges$FibreEnergy[1] & FibreEnergy <= normal_ranges$FibreEnergy[2])),
    (is.na(SodiumPotassiumRatio) | (SodiumPotassiumRatio >= normal_ranges$SodiumPotassiumRatio[1] & SodiumPotassiumRatio <= normal_ranges$SodiumPotassiumRatio[2]))
  ) %>%
  mutate(Hypertension = as.factor(case_when(Systolic >= 120 | Diastolic >= 80 ~ 1,TRUE ~ 0)))

Variable check

variable_data <- data.frame(
  Variable = c("PersonID", "HouseID", "SocialStatus", "Age", "Gender", "BMI", 
               "Systolic", "Diastolic", "SmokeStatus", "AlcoholPercentage", 
               "SodiumPotassiumRatio", "FibreEnergy", "CarbohydrateEnergy", 
               "Identity"),
  Description = c("Selected person identifier", "Household identifier",
                  "ICD10 Conditions - Mini classification",
                  "Age of person", "Sex of person",
                  "Body Mass Index (BMI) - score measured",
                  "Systolic Blood Pressure (mmHg)",
                  "Diastolic Blood Pressure (mmHg)", "Smoking Status",
                  "Average of Percentage of energy from alcohol (Day1) and Percentage of energy from alcohol (Day2)",
                  "The ratio of the average of Sodium (supplement) Day 1 mg and Sodium (supplement) Day 2 mg over the average of Potassium (total) Day 1 mg and Potassium (total) Day 2 mg",
                  "Average of Percentage of energy from fibre (Day1) and Percentage of energy from fibre (Day2)",
                  "Average of Percentage of energy from carbohydrate (Day1) and Percentage of energy from carbohydrate (Day2)",
                  "Indicates whether the individual is part of the indigenous population"))

variable_table <- variable_data %>%
  gt() %>%
  tab_header(title = "All related variables") %>%
  cols_label(Variable = "Variable",Description = "Description") %>%
  tab_caption(caption = md("**Table 4:** Variable Descriptions"))%>%
  fmt_markdown(columns = everything()) %>%
  tab_options(table.width = pct(100))

variable_table

Table 4: Variable Descriptions

All related variables
Variable	Description
PersonID	Selected person identifier
HouseID	Household identifier
SocialStatus	ICD10 Conditions - Mini classification
Age	Age of person
Gender	Sex of person
BMI	Body Mass Index (BMI) - score measured
Systolic	Systolic Blood Pressure (mmHg)
Diastolic	Diastolic Blood Pressure (mmHg)
SmokeStatus	Smoking Status
AlcoholPercentage	Average of Percentage of energy from alcohol (Day1) and Percentage of energy from alcohol (Day2)
SodiumPotassiumRatio	The ratio of the average of Sodium (supplement) Day 1 mg and Sodium (supplement) Day 2 mg over the average of Potassium (total) Day 1 mg and Potassium (total) Day 2 mg
FibreEnergy	Average of Percentage of energy from fibre (Day1) and Percentage of energy from fibre (Day2)
CarbohydrateEnergy	Average of Percentage of energy from carbohydrate (Day1) and Percentage of energy from carbohydrate (Day2)
Identity	Indicates whether the individual is part of the indigenous population

selected_predictors <- selected_data %>%
  dplyr::select(FibreEnergy, CarbohydrateEnergy)

predictor_table <- data.frame(Variable = names(selected_predictors), Type = sapply(selected_predictors, class)) %>%
  gt() %>%
  tab_header(title = "Variable Names and Their Types") %>%
  tab_caption(caption = md("**Table 5**: Predictor variables and their data types"))

predictor_table

Table 5: Predictor variables and their data types

Variable Names and Their Types
Variable	Type
FibreEnergy	numeric
CarbohydrateEnergy	numeric

confounding_vars <- selected_data %>%
  dplyr::select(-FibreEnergy, -CarbohydrateEnergy,-PersonID, -HouseID)

confounding_table <- data.frame(Variable = names(confounding_vars), Type = sapply(confounding_vars, class)) %>%
  gt() %>%
  tab_header(title = "Confounding Variables and Their Types") %>%
  tab_caption(caption = md("**Table 6**: Confounding variables and their data types"))

confounding_table

Table 6: Confounding variables and their data types

Confounding Variables and Their Types
Variable	Type
SocialStatus	factor
Age	integer
Gender	factor
BMI	numeric
Systolic	integer
Diastolic	integer
SmokeStatus	factor
AlcoholPercentage	numeric
SodiumPotassiumRatio	numeric
Identity	factor
Hypertension	factor

format_counts <- function(count, total) {
  paste(count, "(", round(count / total * 100, 1), "%)", sep = "")
}

format_counts_list <- function(values, data, column) {
  sapply(values, function(x) format_counts(sum(data[[column]] == x), nrow(data)))
}

format_mean_sd <- function(mean_val, sd_val) {
  paste(round(mean_val, 2), "(", round(sd_val, 2), ")", sep = "")
}

format_mean_sd_list <- function(data, columns) {
  sapply(columns, function(col) format_mean_sd(mean(data[[col]], na.rm = TRUE), sd(data[[col]], na.rm = TRUE)))
}

get_social_status <- function(data) {
  format_counts_list(1:5, data, "SocialStatus")
}

get_gender <- function(data) {
  format_counts_list(1:2, data, "Gender")
}

get_smoke_status <- function(data) {
  format_counts_list(1:3, data, "SmokeStatus")
}

get_identity <- function(data) {
  format_counts_list(0:1, data, "Identity")
}

get_age <- function(data) {
  c(
    format_counts(sum(data$Age >= 18 & data$Age <= 49), nrow(data)),
    format_counts(sum(data$Age > 50), nrow(data))
  )
}

get_health_metrics <- function(data) {
  format_mean_sd_list(data, c("BMI", "Systolic", "Diastolic", "AlcoholPercentage", "PotassiumSodiumRatio", "FiberEnergy", "CarbonhydrateEnergy"))
}

riskdata <- selected_data %>% filter(Hypertension == 1)
noriskdata <- selected_data %>% filter(Hypertension == 0)

df <- data.frame(
  Variable = c(
    "Lowest 20%, N(%)", "Second quintile, N(%)", "Third quintile, N(%)", "Fourth quintile, N(%)", "Highest 20%*, N(%)",
    "Male, N(%)", "Female, N(%)",
    "Currently Smoker, N(%)", "Ex-Smoker, N(%)", "Never smoked, N(%)",
    "Non Indigenous, N(%)", "Indigenous, N(%)",
    "18-49, N(%)", ">50, N(%)",
    "BMI, kg/m² , mean(SD)", 
    "Systolic, mmHg, mean(SD)", "Diastolic, mmHg, mean(SD)", 
    "Alcohol Percentage, mean(SD)", 
    "Potassium Sodium Ratio, mean(SD)", 
    "Fiber Energy, mean(SD)", 
    "Carbohydrate Energy, mean(SD)"
  ),
  Group = c(
    "SocialStatus", "SocialStatus", "SocialStatus", "SocialStatus", "SocialStatus",
    "Gender", "Gender",
    "SmokeStatus", "SmokeStatus", "SmokeStatus",
    "Identity", "Identity",
    "Age", "Age",
    NA, NA, NA, NA, NA, NA, NA
  ),
  Total = c(
    get_social_status(selected_data),
    get_gender(selected_data),
    get_smoke_status(selected_data),
    get_identity(selected_data),
    get_age(selected_data),
    get_health_metrics(selected_data)
  ),
  HypertensionRisk = c(
    get_social_status(riskdata),
    get_gender(riskdata),
    get_smoke_status(riskdata),
    get_identity(riskdata),
    get_age(riskdata),
    get_health_metrics(riskdata)
  ),
  NoHypertensionRisk = c(
    get_social_status(noriskdata),
    get_gender(noriskdata),
    get_smoke_status(noriskdata),
    get_identity(noriskdata),
    get_age(noriskdata),
    get_health_metrics(noriskdata)
  )
)

total_count <- nrow(selected_data)
risk_count <- nrow(riskdata)
no_risk_count <- nrow(noriskdata)

df %>%
  format_table(ci_brackets = c("(", ")")) %>%
  export_table(format = "html", title = "Table 7: Demographic and Health Characteristics by Hypertension Status")

Table 7: Demographic and Health Characteristics by Hypertension Status
Variable	Total	HypertensionRisk	NoHypertensionRisk
SocialStatus
Lowest 20%, N(%)	2164(27.3%)	1079(28.7%)	1085(26%)
Second quintile, N(%)	1570(19.8%)	769(20.4%)	801(19.2%)
Third quintile, N(%)	1430(18%)	665(17.7%)	765(18.3%)
Fourth quintile, N(%)	1223(15.4%)	572(15.2%)	651(15.6%)
Highest 20%*, N(%)	1545(19.5%)	678(18%)	867(20.8%)
Gender
Male, N(%)	3639(45.9%)	2001(53.2%)	1638(39.3%)
Female, N(%)	4293(54.1%)	1762(46.8%)	2531(60.7%)
SmokeStatus
Currently Smoker, N(%)	2163(27.3%)	1038(27.6%)	1125(27%)
Ex-Smoker, N(%)	2222(28%)	1112(29.6%)	1110(26.6%)
Never smoked, N(%)	3547(44.7%)	1613(42.9%)	1934(46.4%)
Identity
Non Indigenous, N(%)	6290(79.3%)	2947(78.3%)	3343(80.2%)
Indigenous, N(%)	1642(20.7%)	816(21.7%)	826(19.8%)
Age
18-49, N(%)	5174(65.2%)	2149(57.1%)	3025(72.6%)
>50, N(%)	2597(32.7%)	1537(40.8%)	1060(25.4%)
BMI, kg/m² , mean(SD)	27.04(5.28)	28.18(5.33)	25.71(4.9)
Systolic, mmHg, mean(SD)	120.49(15.92)	130.47(13.58)	107.74(7.21)
Diastolic, mmHg, mean(SD)	77.26(9.73)	82.73(8.91)	70.27(5.2)
Alcohol Percentage, mean(SD)	3.3(6.33)	3.84(6.94)	2.81(5.69)
Potassium Sodium Ratio, mean(SD)	NA(NA)	NA(NA)	NA(NA)
Fiber Energy, mean(SD)	NA(NA)	NA(NA)	NA(NA)
Carbohydrate Energy, mean(SD)	NA(NA)	NA(NA)	NA(NA)

Missing Values

selected_data %>% 
  dplyr::select(-PersonID, -HouseID) %>% 
  miss_var_summary() %>% 
  gt() %>%
  tab_header(title = "Missing Value") %>%
  tab_caption(caption = md("Table 8: Missing Values in Selected Data")) %>%
  cols_label(variable = "Variable Name", n_miss = "Missing Count", pct_miss = "Percentage Missing") %>%
  fmt_number(columns = everything(), decimals = 2)

Table 8: Missing Values in Selected Data

Missing Value
Variable Name	Missing Count	Percentage Missing
BMI	1,242.00	15.7
Systolic	1,222.00	15.4
Diastolic	1,222.00	15.4
SocialStatus	0.00	0
Age	0.00	0
Gender	0.00	0
SmokeStatus	0.00	0
AlcoholPercentage	0.00	0
SodiumPotassiumRatio	0.00	0
FibreEnergy	0.00	0
CarbohydrateEnergy	0.00	0
Identity	0.00	0
Hypertension	0.00	0

corr_matrix <- selected_data%>%
  dplyr::select(where(is.numeric)) %>%
  mutate(across(everything(), ~ ifelse(is.na(.), mean(., na.rm = TRUE), .))) %>%
  cor(use = "complete.obs")

ggcorrplot(corr_matrix, lab = TRUE, 
           title = "Figure 2: Correlation Heatmap for ALL", 
           colors = c("blue", "white", "red"),
           tl.cex = 9, lab_size = 4) + theme(legend.position = "none")

Regression Imputation

selected_data$BMI[is.na(selected_data$BMI)] <- lm(BMI ~ Gender + SocialStatus + Age + SodiumPotassiumRatio + FibreEnergy + CarbohydrateEnergy + AlcoholPercentage + SmokeStatus + Identity, data = selected_data, na.action = na.exclude) %>% predict(newdata = selected_data[is.na(selected_data$BMI), ])

selected_data$Systolic[is.na(selected_data$Systolic)] <- lm(Systolic ~ Gender + SocialStatus + Age + BMI + SodiumPotassiumRatio + FibreEnergy + CarbohydrateEnergy + AlcoholPercentage + SmokeStatus + Identity, data = selected_data, na.action = na.exclude) %>% predict(newdata = selected_data[is.na(selected_data$Systolic), ])

selected_data$Diastolic[is.na(selected_data$Diastolic)] <- lm(Diastolic ~ Gender + SocialStatus + Age + BMI + Systolic + SodiumPotassiumRatio + FibreEnergy + CarbohydrateEnergy + AlcoholPercentage + SmokeStatus + Identity, data = selected_data, na.action = na.exclude) %>% predict(newdata = selected_data[is.na(selected_data$Diastolic), ])

vis_miss(selected_data %>% dplyr::select(-PersonID, -HouseID)) +
  ggtitle("Figure 3: Missing Data Visualization of Selected Variables")

Normalization

numeric_data <- selected_data %>% dplyr::select_if(is.numeric)
preprocess_params <- preProcess(numeric_data, method = c("center", "scale"))
scaled_numeric_data <- predict(preprocess_params, numeric_data)
normalized_data <- bind_cols(scaled_numeric_data, selected_data %>% dplyr::select_if(Negate(is.numeric)))

Multicollinearity

selected_data %>%
  dplyr::select(where(is.numeric), -Systolic, -Diastolic) %>%
  ggscatmat() +
  ggtitle("Figure 4: Scatterplot Matrix of Selected Numeric Variables")

Correlation

Scatter plot

fiber_syst_plot <- ggplot(selected_data, aes(x = FibreEnergy, y = Systolic)) +
  geom_point(color = "#FF7F0E", alpha = 0.5) +
  geom_smooth(method = "lm", se = FALSE, color = "#1F77B4", size = 1.5) +
  labs(title = "Figure 5.1: Fiber Energy vs Systolic", x = "Percentage of Energy from Fiber", y = "Systolic(mmHg)") +
  theme_light(base_size = 15) +
  theme(plot.title = element_text(hjust = 0.5))

fiber_dias_plot <- ggplot(selected_data, aes(x = FibreEnergy, y = Diastolic)) +
  geom_point(color = "#FF7F0E", alpha = 0.5) +
  geom_smooth(method = "lm", se = FALSE, color = "#1F77B4", size = 1.5) +
  labs(title = "Figure 5.2: Fiber Energy vs Diastolic", x = "Percentage of Energy from Fiber", y = "Diastolic(mmHg)") +
  theme_light(base_size = 15) +
  theme(plot.title = element_text(hjust = 0.5))

carbo_syst_plot <- ggplot(selected_data, aes(x = CarbohydrateEnergy, y = Systolic)) +
  geom_point(color = "#2CA02C", alpha = 0.5) +
  geom_smooth(method = "lm", se = FALSE, color = "#D62728", size = 1.5) +
  labs(title = "Figure 5.3: Carbohydrate Energy vs Systolic", x = "Percentage of Energy from Carbohydrate", y = "Systolic(mmHg)") +
  theme_light(base_size = 15) +
  theme(plot.title = element_text(hjust = 0.5))

carbo_dias_plot <- ggplot(selected_data, aes(x = CarbohydrateEnergy, y = Diastolic)) +
  geom_point(color = "#2CA02C", alpha = 0.5) +
  geom_smooth(method = "lm", se = FALSE, color = "#D62728", size = 1.5) +
  labs(title = "Figure 5.4: Carbohydrate Energy vs Diastolic", x = "Percentage of Energy from Carbohydrate", y = "Diastolic(mmHg)") +
  theme_light(base_size = 15) +
  theme(plot.title = element_text(hjust = 0.5))

grid.arrange(fiber_syst_plot, fiber_dias_plot, carbo_syst_plot, carbo_dias_plot, nrow = 2)

Linear correlation table

fiber_syst_pearson <- cor(selected_data$FibreEnergy, selected_data$Systolic, method = "pearson")
fiber_syst_spearman <- cor(selected_data$FibreEnergy, selected_data$Systolic, method = "spearman")
fiber_dias_pearson <- cor(selected_data$FibreEnergy, selected_data$Diastolic, method = "pearson")
fiber_dias_spearman <- cor(selected_data$FibreEnergy, selected_data$Diastolic, method = "spearman")

carbo_syst_pearson <- cor(selected_data$CarbohydrateEnergy, selected_data$Systolic, method = "pearson")
carbo_syst_spearman <- cor(selected_data$CarbohydrateEnergy, selected_data$Systolic, method = "spearman")
carbo_dias_pearson <- cor(selected_data$CarbohydrateEnergy, selected_data$Diastolic, method = "pearson")
carbo_dias_spearman <- cor(selected_data$CarbohydrateEnergy, selected_data$Diastolic, method = "spearman")

correlation_table <- data.frame(
  Measure = c("FibreEnergy vs Systolic", "FibreEnergy vs Diastolic",
              "CarbohydrateEnergy vs Systolic", "CarbohydrateEnergy vs Diastolic"),
  Pearson_Correlation = c(fiber_syst_pearson, fiber_dias_pearson, carbo_syst_pearson, carbo_dias_pearson),
  Spearman_Correlation = c(fiber_syst_spearman, fiber_dias_spearman, carbo_syst_spearman, carbo_dias_spearman))

correlation_table %>%
  gt() %>%
  tab_header(
    title = "Table 9: Correlation Between Nutrient Energy Intake and Blood Pressure",
    subtitle = "Correlation values for both Systolic and Diastolic Blood Pressure") %>%
  fmt_number(columns = c(Pearson_Correlation, Spearman_Correlation), decimals = 2)

Table 9: Correlation Between Nutrient Energy Intake and Blood Pressure
Correlation values for both Systolic and Diastolic Blood Pressure
Measure	Pearson_Correlation	Spearman_Correlation
FibreEnergy vs Systolic	0.02	0.02
FibreEnergy vs Diastolic	−0.09	−0.09
CarbohydrateEnergy vs Systolic	−0.04	−0.05
CarbohydrateEnergy vs Diastolic	−0.09	−0.10

Non-linear correlation table

fiber_syst_kendall <- cor(selected_data$FibreEnergy, selected_data$Systolic, method = "kendall")
fiber_syst_mic <- mine(selected_data$FibreEnergy, selected_data$Systolic)$MIC
fiber_dias_kendall <- cor(selected_data$FibreEnergy, selected_data$Diastolic, method = "kendall")
fiber_dias_mic <- mine(selected_data$FibreEnergy, selected_data$Diastolic)$MIC
carbo_syst_kendall <- cor(selected_data$CarbohydrateEnergy, selected_data$Systolic, method = "kendall")
carbo_syst_mic <- mine(selected_data$CarbohydrateEnergy, selected_data$Systolic)$MIC
carbo_dias_kendall <- cor(selected_data$CarbohydrateEnergy, selected_data$Diastolic, method = "kendall")
carbo_dias_mic <- mine(selected_data$CarbohydrateEnergy, selected_data$Diastolic)$MIC

correlation_table_mic_kendall <- data.frame(
  Measure = c("FibreEnergy vs Systolic", "FibreEnergy vs Diastolic",
              "CarbohydrateEnergy vs Systolic", "CarbohydrateEnergy vs Diastolic"),
  Kendall_Correlation = c(fiber_syst_kendall, fiber_dias_kendall, carbo_syst_kendall, carbo_dias_kendall),
  MIC_Correlation = c(fiber_syst_mic, fiber_dias_mic, carbo_syst_mic, carbo_dias_mic))

correlation_table_mic_kendall %>%
  gt() %>%
  tab_header(
    title = "Table 10: Kendall's Tau and MIC Correlation between Nutrient Energy Intake and Blood Pressure") %>%
  fmt_number(columns = c(Kendall_Correlation, MIC_Correlation), decimals = 3)

Table 10: Kendall's Tau and MIC Correlation between Nutrient Energy Intake and Blood Pressure
Measure	Kendall_Correlation	MIC_Correlation
FibreEnergy vs Systolic	0.010	0.044
FibreEnergy vs Diastolic	−0.063	0.058
CarbohydrateEnergy vs Systolic	−0.033	0.059
CarbohydrateEnergy vs Diastolic	−0.070	0.072

Hypothesis testing

median_fiber <- median(selected_data$FibreEnergy)
selected_data$FiberGroup <- ifelse(selected_data$FibreEnergy > median_fiber, "High Fiber", "Low Fiber")

median_carb <- median(selected_data$CarbohydrateEnergy)
selected_data$CarboGroup <- ifelse(selected_data$CarbohydrateEnergy > median_carb, "High Carbo", "Low Carbo")

fiber_systolic_anova <- aov(Systolic ~ FiberGroup, data = selected_data)
fiber_diastolic_anova <- aov(Diastolic ~ FiberGroup, data = selected_data)
carbo_systolic_anova <- aov(Systolic ~ CarboGroup, data = selected_data)
carbo_diastolic_anova <- aov(Diastolic ~ CarboGroup, data = selected_data)

anova_results <- data.frame(
  Measure = c("Fiber vs Systolic", "Fiber vs Diastolic", 
              "Carbohydrate vs Systolic", "Carbohydrate vs Diastolic"),
  F_Statistic = c(summary(fiber_systolic_anova)[[1]][["F value"]][1],
                  summary(fiber_diastolic_anova)[[1]][["F value"]][1],
                  summary(carbo_systolic_anova)[[1]][["F value"]][1],
                  summary(carbo_diastolic_anova)[[1]][["F value"]][1]),
  P_Value = c(summary(fiber_systolic_anova)[[1]][["Pr(>F)"]][1],
              summary(fiber_diastolic_anova)[[1]][["Pr(>F)"]][1],
              summary(carbo_systolic_anova)[[1]][["Pr(>F)"]][1],
              summary(carbo_diastolic_anova)[[1]][["Pr(>F)"]][1])
)

anova_results %>%
  gt() %>%
  tab_header(
    title = "Table 11: ANOVA Test Results: Nutrient Energy Intake vs Blood Pressure"
  ) %>%
  fmt_number(columns = c(F_Statistic, P_Value), decimals = 3) %>%
  cols_label(
    Measure = "Nutrient vs Blood Pressure",
    F_Statistic = "F-Statistic",
    P_Value = "P-Value"
  ) %>%
  tab_options(
    table.font.size = "medium",
    heading.align = "center"
  )

Table 11: ANOVA Test Results: Nutrient Energy Intake vs Blood Pressure
Nutrient vs Blood Pressure	F-Statistic	P-Value
Fiber vs Systolic	0.842	0.359
Fiber vs Diastolic	41.924	0.000
Carbohydrate vs Systolic	5.385	0.020
Carbohydrate vs Diastolic	52.412	0.000

Model comparison

Systolic

set.seed(3888)
k <- 5
repeats <- 3

gam_rmse_systolic <- numeric()
gam_mae_systolic <- numeric()
gam_r2_systolic <- numeric()

mlr_rmse_systolic <- numeric()
mlr_mae_systolic <- numeric()
mlr_r2_systolic <- numeric()

robust_rmse_systolic <- numeric()
robust_mae_systolic <- numeric()
robust_r2_systolic <- numeric()

n <- nrow(normalized_data)

for (r in 1:repeats) {
  folds <- sample(rep(1:k, length.out = n))
  
  for (i in 1:k) {
    train_idx <- which(folds != i)
    test_idx <- which(folds == i)
    
    data_train <- normalized_data[train_idx, ]
    data_test <- normalized_data[test_idx, ]
    
    actual <- data_test$Systolic
    #GAM
    gam_model <- gam(Systolic ~ s(Age) + s(BMI) + SocialStatus + Gender + SmokeStatus +
                       s(AlcoholPercentage) + s(FibreEnergy) + s(CarbohydrateEnergy) + 
                       Identity + s(SodiumPotassiumRatio),
                     data = data_train)
    
    gam_predictions <- predict(gam_model, newdata = data_test)
    
    gam_rmse <- sqrt(mean((gam_predictions - actual)^2))
    gam_mae <- mean(abs(gam_predictions - actual))
    gam_r2 <- cor(gam_predictions, actual)^2
    
    gam_rmse_systolic <- c(gam_rmse_systolic, gam_rmse)
    gam_mae_systolic <- c(gam_mae_systolic, gam_mae)
    gam_r2_systolic <- c(gam_r2_systolic, gam_r2)
    #Multiple Linear Regression
    mlr_model <- lm(Systolic ~ Age + BMI + SocialStatus + Gender + SmokeStatus +
                      AlcoholPercentage + FibreEnergy + CarbohydrateEnergy + 
                      Identity + SodiumPotassiumRatio, data = data_train)
    
    mlr_predictions <- predict(mlr_model, newdata = data_test)
    
    mlr_rmse <- sqrt(mean((mlr_predictions - actual)^2))
    mlr_mae <- mean(abs(mlr_predictions - actual))
    mlr_r2 <- cor(mlr_predictions, actual)^2
    
    mlr_rmse_systolic <- c(mlr_rmse_systolic, mlr_rmse)
    mlr_mae_systolic <- c(mlr_mae_systolic, mlr_mae)
    mlr_r2_systolic <- c(mlr_r2_systolic, mlr_r2)
    #Robust Regression
    robust_model <- rlm(Systolic ~ Age + BMI + SocialStatus + Gender + SmokeStatus +
                          AlcoholPercentage + FibreEnergy + CarbohydrateEnergy + 
                          Identity + SodiumPotassiumRatio, data = data_train)
    
    robust_predictions <- predict(robust_model, newdata = data_test)
    
    robust_rmse <- sqrt(mean((robust_predictions - actual)^2))
    robust_mae <- mean(abs(robust_predictions - actual))
    robust_r2 <- cor(robust_predictions, actual)^2
    
    robust_rmse_systolic <- c(robust_rmse_systolic, robust_rmse)
    robust_mae_systolic <- c(robust_mae_systolic, robust_mae)
    robust_r2_systolic <- c(robust_r2_systolic, robust_r2)
  }
}

mean_systolic_gam_rmse <- mean(gam_rmse_systolic)
mean_systolic_gam_mae <- mean(gam_mae_systolic)
mean_systolic_gam_r2 <- mean(gam_r2_systolic)

mean_systolic_mlr_rmse <- mean(mlr_rmse_systolic)
mean_systolic_mlr_mae <- mean(mlr_mae_systolic)
mean_systolic_mlr_r2 <- mean(mlr_r2_systolic)

mean_systolic_robust_rmse <- mean(robust_rmse_systolic)
mean_systolic_robust_mae <- mean(robust_mae_systolic)
mean_systolic_robust_r2 <- mean(robust_r2_systolic)

Systolic_df <- data.frame(
    Model = c("GAM", "Multiple Linear Regression", "Robust"),
    RMSE = c(mean_systolic_gam_rmse, mean_systolic_mlr_rmse, mean_systolic_robust_rmse),
    MAE = c(mean_systolic_gam_mae, mean_systolic_mlr_mae, mean_systolic_robust_mae),
    R2 = c(mean_systolic_gam_r2, mean_systolic_mlr_r2, mean_systolic_robust_r2)
)

# Define custom colors for the bar plots
colors <- c('#FF5179', '#BF59D9', '#E76BF3')

rmse_range <- c(min(Systolic_df$RMSE) * 0.95, max(Systolic_df$RMSE) * 1.05)
mae_range <- c(min(Systolic_df$MAE) * 0.95, max(Systolic_df$MAE) * 1.05)
r2_range <- c(min(Systolic_df$R2) * 0.95, max(Systolic_df$R2) * 1.05)

Systolic_rmse <- plot_ly(Systolic_df, x = ~Model, y = ~round(RMSE, 3), type = 'bar', color = ~Model, colors = colors,
               text = ~round(RMSE, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'RMSE', range = rmse_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

Systolic_mae <- plot_ly(Systolic_df, x = ~Model, y = ~round(MAE, 3), type = 'bar', color = ~Model, colors = colors,
               text = ~round(MAE, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'MAE', range = mae_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

Systolic_r2 <- plot_ly(Systolic_df, x = ~Model, y = ~round(R2, 3), type = 'bar', color = ~Model, colors = colors,
               text = ~round(R2, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'R²', range = r2_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

SystolicFig <- subplot(Systolic_rmse, Systolic_mae, Systolic_r2, nrows = 1, shareY = FALSE) %>%
  layout(title = list(text = "Figure 6: RMSE [Left], MAE [Middle], R² [Right] for Full Model and Reduced Model(Systolic)", 
                      font = list(size = 12)))

SystolicFig

Diastolic

set.seed(3888)
gam_rmse_diastolic <- numeric()
gam_mae_diastolic <- numeric()
gam_r2_diastolic <- numeric()

mlr_rmse_diastolic <- numeric()
mlr_mae_diastolic <- numeric()
mlr_r2_diastolic <- numeric()

robust_rmse_diastolic <- numeric()
robust_mae_diastolic <- numeric()
robust_r2_diastolic <- numeric()

for (r in 1:repeats) {
  folds <- sample(rep(1:k, length.out = n))
  
  for (i in 1:k) {
    train_idx <- which(folds != i)
    test_idx <- which(folds == i)
    data_train <- normalized_data[train_idx, ]
    data_test <- normalized_data[test_idx, ]
    
    actual <- data_test$Diastolic
    #GAM
    gam_model <- gam(Diastolic ~ s(Age) + s(BMI) + SocialStatus + Gender + SmokeStatus +
                       s(AlcoholPercentage) + s(FibreEnergy) + s(CarbohydrateEnergy) + 
                       Identity + s(SodiumPotassiumRatio),
                     data = data_train)
    
    gam_predictions <- predict(gam_model, newdata = data_test)
    
    gam_rmse <- sqrt(mean((gam_predictions - actual)^2))
    gam_mae <- mean(abs(gam_predictions - actual))
    gam_r2 <- cor(gam_predictions, actual)^2
    
    gam_rmse_diastolic <- c(gam_rmse_diastolic, gam_rmse)
    gam_mae_diastolic <- c(gam_mae_diastolic, gam_mae)
    gam_r2_diastolic <- c(gam_r2_diastolic, gam_r2)
    #Multiple Linear Regression
    mlr_model <- lm(Diastolic ~ Age + BMI + SocialStatus + Gender + SmokeStatus +
                      AlcoholPercentage + FibreEnergy + CarbohydrateEnergy + 
                      Identity + SodiumPotassiumRatio, data = data_train)
    
    mlr_predictions <- predict(mlr_model, newdata = data_test)
    
    mlr_rmse <- sqrt(mean((mlr_predictions - actual)^2))
    mlr_mae <- mean(abs(mlr_predictions - actual))
    mlr_r2 <- cor(mlr_predictions, actual)^2
    
    mlr_rmse_diastolic <- c(mlr_rmse_diastolic, mlr_rmse)
    mlr_mae_diastolic <- c(mlr_mae_diastolic, mlr_mae)
    mlr_r2_diastolic <- c(mlr_r2_diastolic, mlr_r2)
    #Robust Regression
    robust_model <- rlm(Diastolic ~ Age + BMI + SocialStatus + Gender + SmokeStatus +
                          AlcoholPercentage + FibreEnergy + CarbohydrateEnergy + 
                          Identity + SodiumPotassiumRatio, data = data_train)
    
    robust_predictions <- predict(robust_model, newdata = data_test)
    
    robust_rmse <- sqrt(mean((robust_predictions - actual)^2))
    robust_mae <- mean(abs(robust_predictions - actual))
    robust_r2 <- cor(robust_predictions, actual)^2
    
    robust_rmse_diastolic <- c(robust_rmse_diastolic, robust_rmse)
    robust_mae_diastolic <- c(robust_mae_diastolic, robust_mae)
    robust_r2_diastolic <- c(robust_r2_diastolic, robust_r2)
  }
}

mean_diastolic_gam_rmse <- mean(gam_rmse_diastolic)
mean_diastolic_gam_mae <- mean(gam_mae_diastolic)
mean_diastolic_gam_r2 <- mean(gam_r2_diastolic)

mean_diastolic_mlr_rmse <- mean(mlr_rmse_diastolic)
mean_diastolic_mlr_mae <- mean(mlr_mae_diastolic)
mean_diastolic_mlr_r2 <- mean(mlr_r2_diastolic)

mean_diastolic_robust_rmse <- mean(robust_rmse_diastolic)
mean_diastolic_robust_mae <- mean(robust_mae_diastolic)
mean_diastolic_robust_r2 <- mean(robust_r2_diastolic)

Diastolic_df <- data.frame(
    Model = c("GAM", "Multiple Linear Regression", "Robust"),
    RMSE = c(mean_diastolic_gam_rmse, mean_diastolic_mlr_rmse, mean_diastolic_robust_rmse),
    MAE = c(mean_diastolic_gam_mae, mean_diastolic_mlr_mae, mean_diastolic_robust_mae),
    R2 = c(mean_diastolic_gam_r2, mean_diastolic_mlr_r2, mean_diastolic_robust_r2)
)

colors <- c('#FF5179', '#BF59D9', '#E76BF3')

rmse_range <- c(min(Diastolic_df$RMSE) * 0.95, max(Diastolic_df$RMSE) * 1.05)
mae_range <- c(min(Diastolic_df$MAE) * 0.95, max(Diastolic_df$MAE) * 1.05)
r2_range <- c(min(Diastolic_df$R2) * 0.95, max(Diastolic_df$R2) * 1.05)

Diastolic_rmse <- plot_ly(Diastolic_df, x = ~Model, y = ~round(RMSE, 3), type = 'bar', color = ~Model, colors = colors,
               text = ~round(RMSE, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'RMSE', range = rmse_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

Diastolic_mae <- plot_ly(Diastolic_df, x = ~Model, y = ~round(MAE, 3), type = 'bar', color = ~Model, colors = colors,
               text = ~round(MAE, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'MAE', range = mae_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

Diastolic_r2 <- plot_ly(Diastolic_df, x = ~Model, y = ~round(R2, 3), type = 'bar', color = ~Model, colors = colors,
               text = ~round(R2, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'R²', range = r2_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

DiastolicFig <- subplot(
  Diastolic_rmse, 
  Diastolic_mae, 
  Diastolic_r2, 
  nrows = 1, 
  shareY = FALSE
) %>%
  layout(title = list(text = "Figure 7: RMSE [Left], MAE [Middle], R² [Right] for Full Model and Reduced Model(Diastolic)", 
                      font = list(size = 12)))

DiastolicFig

Full model vs Reduced model

set.seed(3888)
k <- 5
repeats <- 3

gam_full_rmse_systolic <- numeric()
gam_full_mae_systolic <- numeric()
gam_full_r2_systolic <- numeric()

gam_reduced_rmse_systolic <- numeric()
gam_reduced_mae_systolic <- numeric()
gam_reduced_r2_systolic <- numeric()

n <- nrow(normalized_data)

for (r in 1:repeats) {
  folds <- sample(rep(1:k, length.out = n))
  
  for (i in 1:k) {
    train_idx <- which(folds != i)
    test_idx <- which(folds == i)
    
    data_train <- normalized_data[train_idx, ]
    data_test <- normalized_data[test_idx, ]
    
    actual <- data_test$Systolic
    #GAM full
    gam_full_model <- gam(Systolic ~ s(Age) + s(BMI) + SocialStatus + Gender + SmokeStatus +
                       s(AlcoholPercentage) + s(FibreEnergy) + s(CarbohydrateEnergy) + 
                       Identity + s(SodiumPotassiumRatio),
                     data = data_train)
    
    gam_full_predictions <- predict(gam_full_model, newdata = data_test)
    
    gam_full_rmse <- sqrt(mean((gam_full_predictions - actual)^2))
    gam_full_mae <- mean(abs(gam_full_predictions - actual))
    gam_full_r2 <- cor(gam_full_predictions, actual)^2
    
    gam_full_rmse_systolic <- c(gam_full_rmse_systolic, gam_full_rmse)
    gam_full_mae_systolic <- c(gam_full_mae_systolic, gam_full_mae)
    gam_full_r2_systolic <- c(gam_full_r2_systolic, gam_full_r2)
    
    #GAM reduced
    gam_reduced_model <- gam(Systolic ~ s(Age) + s(BMI) + SocialStatus + Gender + SmokeStatus +
                       s(AlcoholPercentage) + Identity + s(SodiumPotassiumRatio), data = data_train)
    gam_reduced_predictions <- predict(gam_reduced_model, newdata = data_test)
    
    gam_reduced_rmse <- sqrt(mean((gam_reduced_predictions - actual)^2))
    gam_reduced_mae <- mean(abs(gam_reduced_predictions - actual))
    gam_reduced_r2 <- cor(gam_reduced_predictions, actual)^2
    
    gam_reduced_rmse_systolic <- c(gam_reduced_rmse_systolic, gam_reduced_rmse)
    gam_reduced_mae_systolic <- c(gam_reduced_mae_systolic, gam_reduced_mae)
    gam_reduced_r2_systolic <- c(gam_reduced_r2_systolic, gam_reduced_r2)
  }
}

mean_systolic_gam_full_rmse <- mean(gam_full_rmse_systolic)
mean_systolic_gam_full_mae <- mean(gam_full_mae_systolic)
mean_systolic_gam_full_r2 <- mean(gam_full_r2_systolic)

mean_systolic_gam_reduced_rmse <- mean(gam_reduced_rmse_systolic)
mean_systolic_gam_reduced_mae <- mean(gam_reduced_mae_systolic)
mean_systolic_gam_reduced_r2 <- mean(gam_reduced_r2_systolic)

Systolic_df <- data.frame(
    Model = c("GAM Full", "GAM Reduced"),
    RMSE = c(mean_systolic_gam_full_rmse, mean_systolic_gam_reduced_rmse),
    MAE = c(mean_systolic_gam_full_mae, mean_systolic_gam_reduced_mae),
    R2 = c(mean_systolic_gam_full_r2, mean_systolic_gam_reduced_r2)
)

# Define custom colors for the bar plots
colors <- c('#FF5179', '#BF59D9')

rmse_range <- c(min(Systolic_df$RMSE) * 0.95, max(Systolic_df$RMSE) * 1.05)
mae_range <- c(min(Systolic_df$MAE) * 0.95, max(Systolic_df$MAE) * 1.05)
r2_range <- c(min(Systolic_df$R2) * 0.95, max(Systolic_df$R2) * 1.05)

Systolic_rmse <- plot_ly(Systolic_df, x = ~Model, y = ~round(RMSE, 3), type = 'bar', color = ~Model, colors = colors,
               text = ~round(RMSE, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'RMSE', range = rmse_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

Systolic_mae <- plot_ly(Systolic_df, x = ~Model, y = ~round(MAE, 3), type = 'bar', color = ~Model, colors = colors,
               text = ~round(MAE, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'MAE', range = mae_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

Systolic_r2 <- plot_ly(Systolic_df, x = ~Model, y = ~round(R2, 3), type = 'bar', color = ~Model, colors = colors,
               text = ~round(R2, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'R²', range = r2_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

SystolicFig <- subplot(Systolic_rmse, Systolic_mae, Systolic_r2, nrows = 1, shareY = FALSE, titleX = TRUE) %>%
  layout(title = list(text = "Figure 8: RMSE [Left], MAE [Middle], R² [Right] for Full and Reduced GAM (Systolic)", 
                      font = list(size = 12)))

SystolicFig

set.seed(3888)
k <- 5
repeats <- 3

gam_full_rmse_diastolic <- numeric()
gam_full_mae_diastolic <- numeric()
gam_full_r2_diastolic <- numeric()

gam_reduced_rmse_diastolic <- numeric()
gam_reduced_mae_diastolic <- numeric()
gam_reduced_r2_diastolic <- numeric()

n <- nrow(normalized_data)

for (r in 1:repeats) {
  folds <- sample(rep(1:k, length.out = n))
  
  for (i in 1:k) {
    train_idx <- which(folds != i)
    test_idx <- which(folds == i)
    
    data_train <- normalized_data[train_idx, ]
    data_test <- normalized_data[test_idx, ]
    
    actual <- data_test$Diastolic
    #GAM full
    gam_full_model <- gam(Diastolic ~ s(Age) + s(BMI) + SocialStatus + Gender + SmokeStatus +
                       s(AlcoholPercentage) + s(FibreEnergy) + s(CarbohydrateEnergy) + 
                       Identity + s(SodiumPotassiumRatio),
                     data = data_train)
    
    gam_full_predictions <- predict(gam_full_model, newdata = data_test)
    
    gam_full_rmse <- sqrt(mean((gam_full_predictions - actual)^2))
    gam_full_mae <- mean(abs(gam_full_predictions - actual))
    gam_full_r2 <- cor(gam_full_predictions, actual)^2
    
    gam_full_rmse_diastolic <- c(gam_full_rmse_diastolic, gam_full_rmse)
    gam_full_mae_diastolic <- c(gam_full_mae_diastolic, gam_full_mae)
    gam_full_r2_diastolic <- c(gam_full_r2_diastolic, gam_full_r2)
    
    #GAM reduced
    gam_reduced_model <- gam(Diastolic ~ s(Age) + s(BMI) + SocialStatus + Gender + SmokeStatus +
                       s(AlcoholPercentage) + Identity + s(SodiumPotassiumRatio), data = data_train)
    gam_reduced_predictions <- predict(gam_reduced_model, newdata = data_test)
    
    gam_reduced_rmse <- sqrt(mean((gam_reduced_predictions - actual)^2))
    gam_reduced_mae <- mean(abs(gam_reduced_predictions - actual))
    gam_reduced_r2 <- cor(gam_reduced_predictions, actual)^2
    
    gam_reduced_rmse_diastolic <- c(gam_reduced_rmse_diastolic, gam_reduced_rmse)
    gam_reduced_mae_diastolic <- c(gam_reduced_mae_diastolic, gam_reduced_mae)
    gam_reduced_r2_diastolic <- c(gam_reduced_r2_diastolic, gam_reduced_r2)
  }
}

mean_diastolic_gam_full_rmse <- mean(gam_full_rmse_diastolic)
mean_diastolic_gam_full_mae <- mean(gam_full_mae_diastolic)
mean_diastolic_gam_full_r2 <- mean(gam_full_r2_diastolic)

mean_diastolic_gam_reduced_rmse <- mean(gam_reduced_rmse_diastolic)
mean_diastolic_gam_reduced_mae <- mean(gam_reduced_mae_diastolic)
mean_diastolic_gam_reduced_r2 <- mean(gam_reduced_r2_diastolic)

Diastolic_df <- data.frame(
    Model = c("GAM Full", "GAM Reduced"),
    RMSE = c(mean_diastolic_gam_full_rmse, mean_diastolic_gam_reduced_rmse),
    MAE = c(mean_diastolic_gam_full_mae, mean_diastolic_gam_reduced_mae),
    R2 = c(mean_diastolic_gam_full_r2, mean_diastolic_gam_reduced_r2)
)

# Define custom colors for the bar plots
colors <- c('#FF5179', '#BF59D9')

rmse_range <- c(min(Diastolic_df$RMSE) * 0.95, max(Diastolic_df$RMSE) * 1.05)
mae_range <- c(min(Diastolic_df$MAE) * 0.95, max(Diastolic_df$MAE) * 1.05)
r2_range <- c(min(Diastolic_df$R2) * 0.95, max(Diastolic_df$R2) * 1.05)

Diastolic_rmse <- plot_ly(Diastolic_df, x = ~Model, y = ~round(RMSE, 3), type = 'bar', color = ~Model, colors = colors,
               text = ~round(RMSE, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'RMSE', range = rmse_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

Diastolic_mae <- plot_ly(Diastolic_df, x = ~Model, y = ~round(MAE, 3), type = 'bar', color = ~Model, colors = colors,
               text = ~round(MAE, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'MAE', range = mae_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

Diastolic_r2 <- plot_ly(Diastolic_df, x = ~Model, y = ~round(R2, 3), type = 'bar', color = ~Model, colors = colors,
               text = ~round(R2, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'R²', range = r2_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

DiastolicFig <- subplot(Diastolic_rmse, Diastolic_mae, Diastolic_r2, nrows = 1, shareY = FALSE, titleX = TRUE) %>%
  layout(title = list(text = "Figure 9: RMSE [Left], MAE [Middle], R² [Right] for Full and Reduced GAM (Diastolic)", 
                      font = list(size = 12)))

DiastolicFig

Stratify

Function for CV

k <- 5
repeats <- 3

sys_cv <- function(data) {
  n <- nrow(data)

  folds <- map(1:repeats, ~ sample(rep(1:k, length.out = n)))
  
  results <- map_dfr(folds, function(fold) {
    map_dfr(1:k, function(i) {
      train_idx <- which(fold != i)
      test_idx <- which(fold == i)
      
      data_train <- data[train_idx, ]
      data_test <- data[test_idx, ]
      
      actual <- data_test$Systolic
      
      gam_model <- gam(Systolic ~ s(FibreEnergy) + s(CarbohydrateEnergy), data = data_train)
      
      gam_predictions <- predict(gam_model, newdata = data_test)
      
      tibble(
        rmse = sqrt(mean((gam_predictions - actual)^2)),
        mae = mean(abs(gam_predictions - actual)),
        r2 = cor(gam_predictions, actual)^2
      )
    })
  })

  summary_results <- results %>%
    summarize(
      mean_rmse = mean(rmse),
      mean_mae = mean(mae),
      mean_r2 = mean(r2)
    )
  
  return(summary_results)
}

dia_cv <- function(data) {
  n <- nrow(data)

  folds <- map(1:repeats, ~ sample(rep(1:k, length.out = n)))
  
  results <- map_dfr(folds, function(fold) {
    map_dfr(1:k, function(i) {
      train_idx <- which(fold != i)
      test_idx <- which(fold == i)
      
      data_train <- data[train_idx, ]
      data_test <- data[test_idx, ]
      
      actual <- data_test$Diastolic
      
      gam_model <- gam(Diastolic ~ s(FibreEnergy) + s(CarbohydrateEnergy), data = data_train)
      
      gam_predictions <- predict(gam_model, newdata = data_test)
      
      tibble(
        rmse = sqrt(mean((gam_predictions - actual)^2)),
        mae = mean(abs(gam_predictions - actual)),
        r2 = cor(gam_predictions, actual)^2
      )
    })
  })

  summary_results <- results %>%
    summarize(
      mean_rmse = mean(rmse),
      mean_mae = mean(mae),
      mean_r2 = mean(r2)
    )
  
  return(summary_results)
}

By Age

Group

age1 <- selected_data %>% filter(Age < 35) %>% select_if(is.numeric)

age2 <- selected_data %>%
  filter(Age >= 35 & Age < 65) %>% select_if(is.numeric)

age3 <- selected_data %>% filter(Age >= 65) %>% select_if(is.numeric)

Systolic

set.seed(3888)
results_age1 <- sys_cv(age1)
results_age2 <- sys_cv(age2)
results_age3 <- sys_cv(age3)

gam_rmse_diastolic_ag1 <- results_age1$mean_rmse
gam_mae_diastolic_ag1 <- results_age1$mean_mae
gam_r2_diastolic_ag1 <- results_age1$mean_r2

gam_rmse_diastolic_ag2 <- results_age2$mean_rmse
gam_mae_diastolic_ag2 <- results_age2$mean_mae
gam_r2_diastolic_ag2 <- results_age2$mean_r2

gam_rmse_diastolic_ag3 <- results_age3$mean_rmse
gam_mae_diastolic_ag3 <- results_age3$mean_mae
gam_r2_diastolic_ag3 <- results_age3$mean_r2

Systolic_df <- data.frame(
    Model = c("18-34", "35-64","65-85"),
    RMSE = c(gam_rmse_diastolic_ag1, gam_rmse_diastolic_ag2, gam_rmse_diastolic_ag3),
    MAE = c(gam_mae_diastolic_ag1, gam_mae_diastolic_ag2, gam_mae_diastolic_ag3),
    R2 = c(gam_r2_diastolic_ag1, gam_r2_diastolic_ag2, gam_r2_diastolic_ag3)
)

colors <- c('#FF5179', '#BF59D9', '#E76BF3')

rmse_range <- c(min(Systolic_df$RMSE) * 0.95, max(Systolic_df$RMSE) * 1.05)
mae_range <- c(min(Systolic_df$MAE) * 0.95, max(Systolic_df$MAE) * 1.05)
r2_range <- c(min(Systolic_df$R2) * 0.95, max(Systolic_df$R2) * 1.05)

Systolic_rmse <- plot_ly(Systolic_df, x = ~Model, y = ~round(RMSE, 3), type = 'bar', color = ~Model, colors = colors,
               text = ~round(RMSE, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'RMSE', range = rmse_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

Systolic_mae <- plot_ly(Systolic_df, x = ~Model, y = ~round(MAE, 3), type = 'bar', color = ~Model, colors = colors,
               text = ~round(MAE, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'MAE', range = mae_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

Systolic_r2 <- plot_ly(Systolic_df, x = ~Model, y = ~round(R2, 3), type = 'bar', color = ~Model, colors = colors,
               text = ~round(R2, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'R²', range = r2_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

AgeSystolicFig <- subplot(Systolic_rmse, Systolic_mae, Systolic_r2, nrows = 1, shareY = FALSE) %>%
  layout(title = list(text = "Figure 10: RMSE [Left], MAE [Middle], R² [Right] for Age Groups (Systolic)"), 
                      font = list(size = 12))

AgeSystolicFig

Diastolic

set.seed(3888)
results_age1 <- dia_cv(age1)
results_age2 <- dia_cv(age2)
results_age3 <- dia_cv(age3)

gam_rmse_diastolic_ag1 <- results_age1$mean_rmse
gam_mae_diastolic_ag1 <- results_age1$mean_mae
gam_r2_diastolic_ag1 <- results_age1$mean_r2

gam_rmse_diastolic_ag2 <- results_age2$mean_rmse
gam_mae_diastolic_ag2 <- results_age2$mean_mae
gam_r2_diastolic_ag2 <- results_age2$mean_r2

gam_rmse_diastolic_ag3 <- results_age3$mean_rmse
gam_mae_diastolic_ag3 <- results_age3$mean_mae
gam_r2_diastolic_ag3 <- results_age3$mean_r2

Diastolic_df <- data.frame(
    Model = c("18-34", "35-64","65-85"),
    RMSE = c(gam_rmse_diastolic_ag1, gam_rmse_diastolic_ag2, gam_rmse_diastolic_ag3),
    MAE = c(gam_mae_diastolic_ag1, gam_mae_diastolic_ag2, gam_mae_diastolic_ag3),
    R2 = c(gam_r2_diastolic_ag1, gam_r2_diastolic_ag2, gam_r2_diastolic_ag3)
)

colors <- c('#FF5179', '#BF59D9', '#E76BF3')

rmse_range <- c(min(Diastolic_df$RMSE) * 0.95, max(Diastolic_df$RMSE) * 1.05)
mae_range <- c(min(Diastolic_df$MAE) * 0.95, max(Diastolic_df$MAE) * 1.05)
r2_range <- c(min(Diastolic_df$R2) * 0.95, max(Diastolic_df$R2) * 1.05)

Diastolic_rmse <- plot_ly(Diastolic_df, x = ~Model, y = ~round(RMSE, 3), type = 'bar', color = ~Model, colors = colors,
               text = ~round(RMSE, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'RMSE', range = rmse_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)'
  )

Diastolic_mae <- plot_ly(Diastolic_df, x = ~Model, y = ~round(MAE, 3), type = 'bar', color = ~Model, colors = colors,
               text = ~round(MAE, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'MAE', range = mae_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)'
  )

Diastolic_r2 <- plot_ly(Diastolic_df, x = ~Model, y = ~round(R2, 3), type = 'bar', color = ~Model, colors = colors,
               text = ~round(R2, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'R²', range = r2_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)'
  )

DiastolicFig <- subplot(Diastolic_rmse, Diastolic_mae, Diastolic_r2, nrows = 1, shareY = FALSE) %>%
  layout(title = list(text = "Figure 11: RMSE [Left], MAE [Middle], R² [Right] for Age Groups (Diastolic)"), font = list(size = 12))

DiastolicFig

Partial effect plots - align with the model summary

age_model_1 <- gam(Systolic ~ s(CarbohydrateEnergy) + s(FibreEnergy), data = age1)
age_model_2 <- gam(Systolic ~ s(CarbohydrateEnergy) + s(FibreEnergy), data = age2)
age_model_3 <- gam(Systolic ~ s(CarbohydrateEnergy) + s(FibreEnergy), data = age3)

age_model_4 <- gam(Diastolic ~ s(CarbohydrateEnergy) + s(FibreEnergy), data = age1)
age_model_5 <- gam(Diastolic ~ s(CarbohydrateEnergy) + s(FibreEnergy), data = age2)
age_model_6 <- gam(Diastolic ~ s(CarbohydrateEnergy) + s(FibreEnergy), data = age3)

FibreEnergy

par(mfrow = c(3, 2))

fiber_color <- "blue"

plot(age_model_1, select = 2, seWithMean = TRUE, shade = TRUE, col = fiber_color, main = "Figure 12.1 18-34 Group - Systolic - FibreEnergy")
plot(age_model_4, select = 2, seWithMean = TRUE, shade = TRUE, col = fiber_color, main = "Figure 12.4 18-34 Group - Diastolic - FibreEnergy")
plot(age_model_2, select = 2, seWithMean = TRUE, shade = TRUE, col = fiber_color, main = "Figure 12.2 35-64 Group - Systolic - FibreEnergy")
plot(age_model_5, select = 2, seWithMean = TRUE, shade = TRUE, col = fiber_color, main = "Figure 12.5 35-64 Group - Diastolic - FibreEnergy")
plot(age_model_3, select = 2, seWithMean = TRUE, shade = TRUE, col = fiber_color, main = "Figure 12.3 65-85 Group - Systolic - FibreEnergy")
plot(age_model_6, select = 2, seWithMean = TRUE, shade = TRUE, col = fiber_color, main = "Figure 12.6 65-85 Group - Diastolic - FibreEnergy")

par(mfrow = c(1, 1))

CarbohydrateEnergy

par(mfrow = c(3, 2))
carb_color <- "red"

plot(age_model_1, select = 1, seWithMean = TRUE, shade = TRUE, col = carb_color, main = "Figure 13.1 18-34 Group - Systolic - CarbohydrateEnergy")
plot(age_model_4, select = 1, seWithMean = TRUE, shade = TRUE, col = carb_color, main = "Figure 13.2 18-34 Group - Diastolic - CarbohydrateEnergy")
plot(age_model_2, select = 1, seWithMean = TRUE, shade = TRUE, col = carb_color, main = "Figure 13.3 35-64 Group - Systolic - CarbohydrateEnergy")
plot(age_model_5, select = 1, seWithMean = TRUE, shade = TRUE, col = carb_color, main = "Figure 13.4 35-64 Group - Diastolic - CarbohydrateEnergy")
plot(age_model_3, select = 1, seWithMean = TRUE, shade = TRUE, col = carb_color, main = "Figure 13.5 65-85 Group - Systolic - CarbohydrateEnergy")
plot(age_model_6, select = 1, seWithMean = TRUE, shade = TRUE, col = carb_color, main = "Figure 13.6 65-85 Group - Diastolic - CarbohydrateEnergy")

par(mfrow = c(1, 1))

By BMI

Group

light <- selected_data %>% filter(BMI < 24.9)

normal <- selected_data %>% filter(BMI >= 24.9 & BMI <= 29.9)

heavy <- selected_data %>% filter(BMI > 29.9)

Systolic

set.seed(3888)
results_light <- sys_cv(light)
results_normal <- sys_cv(normal)
results_heavy <- sys_cv(heavy)

gam_rmse_systolic_light <- results_light$mean_rmse
gam_mae_systolic_light <- results_light$mean_mae
gam_r2_systolic_light <- results_light$mean_r2

gam_rmse_systolic_normal <- results_normal$mean_rmse
gam_mae_systolic_normal <- results_normal$mean_mae
gam_r2_systolic_normal <- results_normal$mean_r2

gam_rmse_systolic_heavy <- results_heavy$mean_rmse
gam_mae_systolic_heavy <- results_heavy$mean_mae
gam_r2_systolic_heavy <- results_heavy$mean_r2

Systolic_df <- data.frame(
    Group = c("Light BMI", "Normal BMI", "Heavy BMI"),
    RMSE = c(gam_rmse_systolic_light, gam_rmse_systolic_normal, gam_rmse_systolic_heavy),
    MAE = c(gam_mae_systolic_light, gam_mae_systolic_normal, gam_mae_systolic_heavy),
    R2 = c(gam_r2_systolic_light, gam_r2_systolic_normal, gam_r2_systolic_heavy)
)

colors <- c('#FF5179', '#BF59D9', '#E76BF3')

rmse_range <- c(min(Systolic_df$RMSE) * 0.95, max(Systolic_df$RMSE) * 1.05)
mae_range <- c(min(Systolic_df$MAE) * 0.95, max(Systolic_df$MAE) * 1.05)
r2_range <- c(min(Systolic_df$R2) * 0.95, max(Systolic_df$R2) * 1.05)

Systolic_rmse <- plot_ly(Systolic_df, x = ~Group, y = ~round(RMSE, 3), type = 'bar', color = ~Group, colors = colors,
               text = ~round(RMSE, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'RMSE', range = rmse_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

Systolic_mae <- plot_ly(Systolic_df, x = ~Group, y = ~round(MAE, 3), type = 'bar', color = ~Group, colors = colors,
               text = ~round(MAE, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'MAE', range = mae_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

Systolic_r2 <- plot_ly(Systolic_df, x = ~Group, y = ~round(R2, 3), type = 'bar', color = ~Group, colors = colors,
               text = ~round(R2, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'R²', range = r2_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

SystolicFig <- subplot(Systolic_rmse, Systolic_mae, Systolic_r2, nrows = 1, shareY = FALSE) %>%
  layout(title = list(text = "Figure 14: RMSE [Left], MAE [Middle], R² [Right] for BMI Groups (Systolic)"), font = list(size = 12))

SystolicFig

Diastolic

set.seed(3888)
results_light_dias <- dia_cv(light)
results_normal_dias <- dia_cv(normal)
results_heavy_dias <- dia_cv(heavy)

gam_rmse_diastolic_light <- results_light_dias$mean_rmse
gam_mae_diastolic_light <- results_light_dias$mean_mae
gam_r2_diastolic_light <- results_light_dias$mean_r2

gam_rmse_diastolic_normal <- results_normal_dias$mean_rmse
gam_mae_diastolic_normal <- results_normal_dias$mean_mae
gam_r2_diastolic_normal <- results_normal_dias$mean_r2

gam_rmse_diastolic_heavy <- results_heavy_dias$mean_rmse
gam_mae_diastolic_heavy <- results_heavy_dias$mean_mae
gam_r2_diastolic_heavy <- results_heavy_dias$mean_r2

Diastolic_df <- data.frame(
    Group = c("Light BMI", "Normal BMI", "Heavy BMI"),
    RMSE = c(gam_rmse_diastolic_light, gam_rmse_diastolic_normal, gam_rmse_diastolic_heavy),
    MAE = c(gam_mae_diastolic_light, gam_mae_diastolic_normal, gam_mae_diastolic_heavy),
    R2 = c(gam_r2_diastolic_light, gam_r2_diastolic_normal, gam_r2_diastolic_heavy)
)

colors <- c('#FF5179', '#BF59D9', '#E76BF3')

rmse_range <- c(min(Diastolic_df$RMSE) * 0.95, max(Diastolic_df$RMSE) * 1.05)
mae_range <- c(min(Diastolic_df$MAE) * 0.95, max(Diastolic_df$MAE) * 1.05)
r2_range <- c(min(Diastolic_df$R2) * 0.95, max(Diastolic_df$R2) * 1.05)

Diastolic_rmse <- plot_ly(Diastolic_df, x = ~Group, y = ~round(RMSE, 3), type = 'bar', color = ~Group, colors = colors,
               text = ~round(RMSE, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'RMSE', range = rmse_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

Diastolic_mae <- plot_ly(Diastolic_df, x = ~Group, y = ~round(MAE, 3), type = 'bar', color = ~Group, colors = colors,
               text = ~round(MAE, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'MAE', range = mae_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

Diastolic_r2 <- plot_ly(Diastolic_df, x = ~Group, y = ~round(R2, 3), type = 'bar', color = ~Group, colors = colors,
               text = ~round(R2, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'R²', range = r2_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

DiastolicFig <- subplot(Diastolic_rmse, Diastolic_mae, Diastolic_r2, nrows = 1, shareY = FALSE) %>%
  layout(title = list(text = "Figure 15: RMSE [Left], MAE [Middle], R² [Right] for BMI Groups (Diastolic)"), font = list(size = 12))

DiastolicFig

Partial effect plots - align with the model summary

bmi_model_1 <- gam(Systolic ~ s(CarbohydrateEnergy) + s(FibreEnergy), data = light)
bmi_model_2 <- gam(Systolic ~ s(CarbohydrateEnergy) + s(FibreEnergy), data = normal)
bmi_model_3 <- gam(Systolic ~ s(CarbohydrateEnergy) + s(FibreEnergy), data = heavy)

bmi_model_4 <- gam(Diastolic ~ s(CarbohydrateEnergy) + s(FibreEnergy), data = light)
bmi_model_5 <- gam(Diastolic ~ s(CarbohydrateEnergy) + s(FibreEnergy), data = normal)
bmi_model_6 <- gam(Diastolic ~ s(CarbohydrateEnergy) + s(FibreEnergy), data = heavy)

FibreEnergy

par(mfrow = c(3, 2))

fiber_color <- "blue"

plot(bmi_model_1, select = 2, seWithMean = TRUE, shade = TRUE, col = fiber_color, main = "Figure 16.1 Low BMI Group - Systolic - FibreEnergy")
plot(bmi_model_4, select = 2, seWithMean = TRUE, shade = TRUE, col = fiber_color, main = "Figure 16.4 Low BMI Group - Diastolic - FibreEnergy")
plot(bmi_model_2, select = 2, seWithMean = TRUE, shade = TRUE, col = fiber_color, main = "Figure 16.2 Normal BMI Group - Systolic - FibreEnergy")
plot(bmi_model_5, select = 2, seWithMean = TRUE, shade = TRUE, col = fiber_color, main = "Figure 16.5 Normal BMI Group - Diastolic - FibreEnergy")
plot(bmi_model_3, select = 2, seWithMean = TRUE, shade = TRUE, col = fiber_color, main = "Figure 16.3 High BMI Group - Systolic - FibreEnergy")
plot(bmi_model_6, select = 2, seWithMean = TRUE, shade = TRUE, col = fiber_color, main = "Figure 16.6 High BMI Group - Diastolic - FibreEnergy")

par(mfrow = c(1, 1))

CarbohydrateEnergy

par(mfrow = c(3, 2))
carb_color <- "red"

plot(bmi_model_1, select = 1, seWithMean = TRUE, shade = TRUE, col = carb_color, main = "Figure 17.1 Low BMI Group - Systolic - CarbohydrateEnergy")
plot(bmi_model_4, select = 1, seWithMean = TRUE, shade = TRUE, col = carb_color, main = "Figure 17.2 Low BMI Group - Diastolic - CarbohydrateEnergy")
plot(bmi_model_2, select = 1, seWithMean = TRUE, shade = TRUE, col = carb_color, main = "Figure 17.3 Normal BMI Group - Systolic - CarbohydrateEnergy")
plot(bmi_model_5, select = 1, seWithMean = TRUE, shade = TRUE, col = carb_color, main = "Figure 17.4 Normal BMI Group - Diastolic - CarbohydrateEnergy")
plot(bmi_model_3, select = 1, seWithMean = TRUE, shade = TRUE, col = carb_color, main = "Figure 17.5 High BMI Group - Systolic - CarbohydrateEnergy")
plot(bmi_model_6, select = 1, seWithMean = TRUE, shade = TRUE, col = carb_color, main = "Figure 17.6 High BMI Group - Diastolic - CarbohydrateEnergy")

par(mfrow = c(1, 1))

By AlcoholPercentage

Group

nonalcohol <- selected_data %>% filter(AlcoholPercentage == 0)

alcohol <- selected_data %>% filter(AlcoholPercentage > 0)

Systolic

set.seed(3888)
results_nonalcohol <- sys_cv(nonalcohol)
results_alcohol <- sys_cv(alcohol)

gam_rmse_systolic_nonalcohol <- results_nonalcohol$mean_rmse
gam_mae_systolic_nonalcohol <- results_nonalcohol$mean_mae
gam_r2_systolic_nonalcohol <- results_nonalcohol$mean_r2

gam_rmse_systolic_alcohol <- results_alcohol$mean_rmse
gam_mae_systolic_alcohol <- results_alcohol$mean_mae
gam_r2_systolic_alcohol <- results_alcohol$mean_r2

Systolic_alcohol_df <- data.frame(
    Group = c("Non-Alcohol", "Alcohol"),
    RMSE = c(gam_rmse_systolic_nonalcohol, gam_rmse_systolic_alcohol),
    MAE = c(gam_mae_systolic_nonalcohol, gam_mae_systolic_alcohol),
    R2 = c(gam_r2_systolic_nonalcohol, gam_r2_systolic_alcohol)
)

colors_alcohol <- c('#FF5179', '#BF59D9')

rmse_range_alcohol <- c(min(Systolic_alcohol_df$RMSE) * 0.95, max(Systolic_alcohol_df$RMSE) * 1.05)
mae_range_alcohol <- c(min(Systolic_alcohol_df$MAE) * 0.95, max(Systolic_alcohol_df$MAE) * 1.05)
r2_range_alcohol <- c(min(Systolic_alcohol_df$R2) * 0.95, max(Systolic_alcohol_df$R2) * 1.05)

Systolic_rmse_alcohol <- plot_ly(Systolic_alcohol_df, x = ~Group, y = ~round(RMSE, 3), type = 'bar', color = ~Group, colors = colors_alcohol,
               text = ~round(RMSE, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'RMSE', range = rmse_range_alcohol),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

Systolic_mae_alcohol <- plot_ly(Systolic_alcohol_df, x = ~Group, y = ~round(MAE, 3), type = 'bar', color = ~Group, colors = colors_alcohol,
               text = ~round(MAE, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'MAE', range = mae_range_alcohol),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

Systolic_r2_alcohol <- plot_ly(Systolic_alcohol_df, x = ~Group, y = ~round(R2, 3), type = 'bar', color = ~Group, colors = colors_alcohol,
               text = ~round(R2, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'R²', range = r2_range_alcohol),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

SystolicFig_alcohol <- subplot(Systolic_rmse_alcohol, Systolic_mae_alcohol, Systolic_r2_alcohol, nrows = 1, shareY = FALSE) %>%
  layout(title = list(text = "Figure 18: RMSE [Left], MAE [Middle], R² [Right] for Alcohol Groups (Systolic)"), 
                      font = list(size = 12))

SystolicFig_alcohol

Diastolic

set.seed(3888)
results_nonalcohol_dias <- dia_cv(nonalcohol)
results_alcohol_dias <- dia_cv(alcohol)

gam_rmse_diastolic_nonalcohol <- results_nonalcohol_dias$mean_rmse
gam_mae_diastolic_nonalcohol <- results_nonalcohol_dias$mean_mae
gam_r2_diastolic_nonalcohol <- results_nonalcohol_dias$mean_r2

gam_rmse_diastolic_alcohol <- results_alcohol_dias$mean_rmse
gam_mae_diastolic_alcohol <- results_alcohol_dias$mean_mae
gam_r2_diastolic_alcohol <- results_alcohol_dias$mean_r2

Diastolic_alcohol_df <- data.frame(
    Group = c("Non-Alcohol", "Alcohol"),
    RMSE = c(gam_rmse_diastolic_nonalcohol, gam_rmse_diastolic_alcohol),
    MAE = c(gam_mae_diastolic_nonalcohol, gam_mae_diastolic_alcohol),
    R2 = c(gam_r2_diastolic_nonalcohol, gam_r2_diastolic_alcohol)
)

colors_alcohol <- c('#FF5179', '#BF59D9')

rmse_range_alcohol <- c(min(Diastolic_alcohol_df$RMSE) * 0.95, max(Diastolic_alcohol_df$RMSE) * 1.05)
mae_range_alcohol <- c(min(Diastolic_alcohol_df$MAE) * 0.95, max(Diastolic_alcohol_df$MAE) * 1.05)
r2_range_alcohol <- c(min(Diastolic_alcohol_df$R2) * 0.95, max(Diastolic_alcohol_df$R2) * 1.05)

Diastolic_rmse_alcohol <- plot_ly(Diastolic_alcohol_df, x = ~Group, y = ~round(RMSE, 3), type = 'bar', color = ~Group, colors = colors_alcohol,
               text = ~round(RMSE, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'RMSE', range = rmse_range_alcohol),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

Diastolic_mae_alcohol <- plot_ly(Diastolic_alcohol_df, x = ~Group, y = ~round(MAE, 3), type = 'bar', color = ~Group, colors = colors_alcohol,
               text = ~round(MAE, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'MAE', range = mae_range_alcohol),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

Diastolic_r2_alcohol <- plot_ly(Diastolic_alcohol_df, x = ~Group, y = ~round(R2, 3), type = 'bar', color = ~Group, colors = colors_alcohol,
               text = ~round(R2, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'R²', range = r2_range_alcohol),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

DiastolicFig_alcohol <- subplot(Diastolic_rmse_alcohol, Diastolic_mae_alcohol, Diastolic_r2_alcohol, nrows = 1, shareY = FALSE) %>%
  layout(title = list(text = "Figure 19: RMSE [Left], MAE [Middle], R² [Right] for Alcohol Groups (Diastolic)"), 
                      font = list(size = 12))

DiastolicFig_alcohol

Partial effect plots - align with the model summary

alc_model_1 <- gam(Systolic ~ s(CarbohydrateEnergy) + s(FibreEnergy), data = nonalcohol)
alc_model_2 <- gam(Systolic ~ s(CarbohydrateEnergy) + s(FibreEnergy), data = alcohol)
alc_model_3 <- gam(Diastolic ~ s(CarbohydrateEnergy) + s(FibreEnergy), data = nonalcohol)
alc_model_4 <- gam(Diastolic ~ s(CarbohydrateEnergy) + s(FibreEnergy), data = alcohol)

FibreEnergy

par(mfrow = c(2, 2))

fiber_color <- "blue"

plot(alc_model_1, select = 2, seWithMean = TRUE, shade = TRUE, col = fiber_color, main = "Figure 20.1 non-alcohol Group - Systolic - FibreEnergy")
plot(alc_model_3, select = 2, seWithMean = TRUE, shade = TRUE, col = fiber_color, main = "Figure 20.4 non-alcohol Group - Diastolic - FibreEnergy")
plot(alc_model_2, select = 2, seWithMean = TRUE, shade = TRUE, col = fiber_color, main = "Figure 20.2 alcohol Group - Systolic - FibreEnergy")
plot(alc_model_4, select = 2, seWithMean = TRUE, shade = TRUE, col = fiber_color, main = "Figure 20.5 alcohol Group - Diastolic - FibreEnergy")

par(mfrow = c(1, 1))

CarbohydrateEnergy

par(mfrow = c(2, 2))
carb_color <- "red"

plot(alc_model_1, select = 1, seWithMean = TRUE, shade = TRUE, col = carb_color, main = "Figure 21.1 non-alcohol Group - Systolic - CarbohydrateEnergy")
plot(alc_model_3, select = 1, seWithMean = TRUE, shade = TRUE, col = carb_color, main = "Figure 21.2 non-alcohol Group - Diastolic - CarbohydrateEnergy")
plot(alc_model_2, select = 1, seWithMean = TRUE, shade = TRUE, col = carb_color, main = "Figure 21.3 alcohol Group - Systolic - CarbohydrateEnergy")
plot(alc_model_4, select = 1, seWithMean = TRUE, shade = TRUE, col = carb_color, main = "Figure 21.4 alcohol Group - Diastolic - CarbohydrateEnergy")

par(mfrow = c(1, 1))

By Gender

Group

male <- selected_data %>% filter(Gender == 1) %>% select_if(is.numeric)

female <- selected_data %>% filter(Gender == 2) %>% select_if(is.numeric)

Systolic

set.seed(3888)
results_male <- sys_cv(male)
results_female <- sys_cv(female)

gam_rmse_diastolic_male <- results_male$mean_rmse
gam_mae_diastolic_male <- results_male$mean_mae
gam_r2_diastolic_male <- results_male$mean_r2

gam_rmse_diastolic_female <- results_female$mean_rmse
gam_mae_diastolic_female <- results_female$mean_mae
gam_r2_diastolic_female <- results_female$mean_r2

Systolic_gender_df <- data.frame(
    Model = c("Male", "Female"),
    RMSE = c(gam_rmse_diastolic_male, gam_rmse_diastolic_female),
    MAE = c(gam_mae_diastolic_male, gam_mae_diastolic_female),
    R2 = c(gam_r2_diastolic_male, gam_r2_diastolic_female)
)

colors <- c('#FF5179', '#BF59D9')

rmse_range <- c(min(Systolic_gender_df$RMSE) * 0.95, max(Systolic_gender_df$RMSE) * 1.05)
mae_range <- c(min(Systolic_gender_df$MAE) * 0.95, max(Systolic_gender_df$MAE) * 1.05)
r2_range <- c(min(Systolic_gender_df$R2) * 0.95, max(Systolic_gender_df$R2) * 1.05)

Systolic_rmse <- plot_ly(Systolic_gender_df, x = ~Model, y = ~round(RMSE, 3), type = 'bar', color = ~Model, colors = colors,
               text = ~round(RMSE, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'RMSE', range = rmse_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

Systolic_mae <- plot_ly(Systolic_gender_df, x = ~Model, y = ~round(MAE, 3), type = 'bar', color = ~Model, colors = colors,
               text = ~round(MAE, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'MAE', range = mae_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

Systolic_r2 <- plot_ly(Systolic_gender_df, x = ~Model, y = ~round(R2, 3), type = 'bar', color = ~Model, colors = colors,
               text = ~round(R2, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'R²', range = r2_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)')

GenderSystolicFig <- subplot(Systolic_rmse, Systolic_mae, Systolic_r2, nrows = 1, shareY = FALSE) %>%
  layout(title = list(text = "Figure 22: RMSE [Left], MAE [Middle], R² [Right] for Gender Groups (Systolic)"), 
                      font = list(size = 12))

GenderSystolicFig

Diastolic

set.seed(3888)

results_male <- dia_cv(male)
results_female <- dia_cv(female)

gam_rmse_diastolic_male <- results_male$mean_rmse
gam_mae_diastolic_male <- results_male$mean_mae
gam_r2_diastolic_male <- results_male$mean_r2

gam_rmse_diastolic_female <- results_female$mean_rmse
gam_mae_diastolic_female <- results_female$mean_mae
gam_r2_diastolic_female <- results_female$mean_r2

Diastolic_gender_df <- data.frame(
    Model = c("Male", "Female"),
    RMSE = c(gam_rmse_diastolic_male, gam_rmse_diastolic_female),
    MAE = c(gam_mae_diastolic_male, gam_mae_diastolic_female),
    R2 = c(gam_r2_diastolic_male, gam_r2_diastolic_female)
)

colors <- c('#FF5179', '#BF59D9')

rmse_range <- c(min(Diastolic_gender_df$RMSE) * 0.95, max(Diastolic_gender_df$RMSE) * 1.05)
mae_range <- c(min(Diastolic_gender_df$MAE) * 0.95, max(Diastolic_gender_df$MAE) * 1.05)
r2_range <- c(min(Diastolic_gender_df$R2) * 0.95, max(Diastolic_gender_df$R2) * 1.05)

Diastolic_rmse <- plot_ly(Diastolic_gender_df, x = ~Model, y = ~round(RMSE, 3), type = 'bar', color = ~Model, colors = colors,
               text = ~round(RMSE, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'RMSE', range = rmse_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)'
  )

Diastolic_mae <- plot_ly(Diastolic_gender_df, x = ~Model, y = ~round(MAE, 3), type = 'bar', color = ~Model, colors = colors,
               text = ~round(MAE, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'MAE', range = mae_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)'
  )

Diastolic_r2 <- plot_ly(Diastolic_gender_df, x = ~Model, y = ~round(R2, 3), type = 'bar', color = ~Model, colors = colors,
               text = ~round(R2, 3), textposition = 'outside', textfont = list(color = 'black')) %>%
  layout(
    yaxis = list(title = 'R²', range = r2_range),
    showlegend = FALSE,
    paper_bgcolor = 'rgba(245, 246, 249, 1)',
    plot_bgcolor = 'rgba(0, 0, 0, 0)'
  )

DiastolicFig <- subplot(Diastolic_rmse, Diastolic_mae, Diastolic_r2, nrows = 1, shareY = FALSE) %>%
  layout(title = list(text = "Figure 23: RMSE [Left], MAE [Middle], R² [Right] for Gender Groups (Diastolic)"), 
                      font = list(size = 12))

DiastolicFig

Partial effect plots - align with the model summary

systolic_male <- gam(Systolic ~ s(CarbohydrateEnergy) + s(FibreEnergy), data = male)
systolic_female <- gam(Systolic ~ s(CarbohydrateEnergy) + s(FibreEnergy), data = female)

diastolic_male <- gam(Diastolic ~ s(CarbohydrateEnergy) + s(FibreEnergy), data = male)
diastolic_female <- gam(Diastolic ~ s(CarbohydrateEnergy) + s(FibreEnergy), data = female)

FibreEnergy

par(mfrow = c(2, 2))

fiber_color <- "blue"

plot(systolic_male, select = 2, seWithMean = TRUE, shade = TRUE, col = fiber_color, main = "Figure 24.1 systolic_male Group - Systolic - FibreEnergy")
plot(diastolic_male, select = 2, seWithMean = TRUE, shade = TRUE, col = fiber_color, main = "Figure 24.2  diastolic_male Group - Diastolic - FibreEnergy")
plot(systolic_female, select = 2, seWithMean = TRUE, shade = TRUE, col = fiber_color, main = "Figure 24.3 systolic_female Group - Systolic - FibreEnergy")
plot(diastolic_female, select = 2, seWithMean = TRUE, shade = TRUE, col = fiber_color, main = "Figure 24.4 distolic_female Group - Diastolic - FibreEnergy")

par(mfrow = c(1, 1))

CarbohydrateEnergy

par(mfrow = c(2, 2))
carb_color <- "red"

plot(systolic_male, select = 1, seWithMean = TRUE, shade = TRUE, col = carb_color, main = "Figure 25.1 systolic_male Group - Systolic - CarbohydrateEnergy")
plot(diastolic_male, select = 1, seWithMean = TRUE, shade = TRUE, col = carb_color, main = "Figure 25.2 diastolic_male Group - Systolic - CarbohydrateEnergy")
plot(systolic_female, select = 1, seWithMean = TRUE, shade = TRUE, col = carb_color, main = "Figure 25.3 systolic_female Group - Diastolic - CarbohydrateEnergy")
plot(diastolic_female, select = 1, seWithMean = TRUE, shade = TRUE, col = carb_color, main = "Figure 25.4 diastolic_female Group - Diastolic - CarbohydrateEnergy")

par(mfrow = c(1, 1))